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ABSTRACT 


By using currently accepted constants of geodesy, notably 
several different sets of oblateness coefficients, seven models of the 
geometry of the earth and its gravitational field, including the cur- 
rently valid NASA Standard Model, are numerically derived and compared 
with each other and with a newly developed geoid based on recent meas- 
urements of the secular and periodic perturbations of artificial 
satellite orbits. 
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DEFINITION OF SYMBOLS 


Symbol 
As, a 


Ba 

^n,m 

f 

£e 

8e 

8e 


GM 




Jn 
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n,m 


J 

H 

D 
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Definition 

coefficients in the series expansion for the geoid radius 

coefficient of the sectorial harmonic in the geoid 
expans ion 

coefficient in the series expansion of the geoid radius 
numerical coefficients in the geopotential function 
flattening or oblateness 

ellipticity (or flattening) of the equator 
gravity at the equator 
mean equatorial gravity 

gravitational parameter of the earth (= (a) 
coefficients of zonal harmonics in geopotential 
coefficients of tesseral harmonics in geopotential 
coefficients of sectorial harmonics in geopotential 

oblateness coefficients (older notation) 

centrifugal factor at the equator 
also; index in series notation 


P™(sin 0 ) associated Legendre polynomials of degree n and order m 
n 


R 

R 

R. 


(sin 0 ) 


Legendre polynomials of nth degree 
equatorial radius of the earth 
mean equatorial radius of the earth 
polar radius of the earth 


V 



Symbol 

^pN 

%S 

^eo 

Re, 

r 

U 

A 


a 

0 

X 

CaJ 


DEFINITION OF SYMBOLS (Continued) 

Definition 

north polar radius of the earth 

south polar radius of the earth 
longest equatorial radius of the earth 
shortest equatorial radius of the earth 
radial distance from the center of the earth 
potential function 

geographic longitude, positive east of Greenwich 

longitude of the longest equatorial radius (semi-major 
axis) 

gravitational parameter of the earth 

factor in equatorial gravity acceleration equation 

geocentric latitude 

factor defined by equation (22) 

angular velocity of the earth's rotation 

centrifugal factor at the equator 

perturbation parameter in potential function 
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SURVEY AND COMPARATIVE ANALYSIS OF CURRENT GEOPHYSICAL MODELS 


SUMMARY 

By using currently accepted constants of geodesy, especially several 
different sets of oblateness coefficients, seven models of the geometry of 
the earth and its gravitational field are numerically derived and compared 
with each other and with a newly developed geopotential based on recent 
measurements of the secular and periodic perturbations of artificial 
satellite orbits. 

Conclusions are drawn which show, among other things, that the cur- 
rently used NASA Standard Model deviates significantly from the "best" 
geoid, and that there are better analytical, as well as numerical, 
models by which Improved approximation can be achieved. 

Also, it is shown that improvements in the numerical values of the 
used oblateness coefficients, rather than in their number, appear to 
have little influence on the values of the gravity, but that they affect 
the geometries significantly. Increasing the number of zonal harmonics 
coeff icients ,, however , can lead to considerable improvements in the 
gravity model, at least in the range studied (up to and including Jx 4 )» 


I. INTRODUCTION 


In the past a number of theories of varying complexity have been 
developed to describe the figure and the gravitation field of the earth 
for use in geodetic, astrodynamic, and astronomical calculations. 

In scrutinizing these theories, a number of discrepancies are mani- 
fest which appear to be caused by two factors; i.e., (1) inherent incon- 
sistencies in the development of some of the theories, and (2) differences 
in gravitational and geometrical parameters used by their authors. 

With regard to the first consideration, while the earth's gravita- 
tional potential (i.e., the potential energy of the earth in reflation to 
the position relative to the earth) is commonly expressed by the poten- 
tial function as a series of spherical harmonics, as adopted by the 
International Astronomical Union (lAU)^*, the shape of the earth 


Superscripts 


in -text indicate references at end of report. 



introduced in the geopotential at the earth's surface with the simul- 
taneous assumption that the potential along the surface of the geoid 
is constant, is that of an ellipsoid which approximates the actual 
geoid. The flattening of the ellipsoid is usually determined approxi- 
matively by the astrogeodetic methods of geometrical geodesy (sometimes 
also by the less accurate techniques of gravimetry), using variants of 
Clairaut's theorem. However, even if higher-order extensions of the 
latter are used, they remain only approximations, which become increas- 
ingly inaccurate for increasing deviation of the geopotential from that 
of a regular ellipsoid. 

The second consideration refers to the fact that, with our increas- 
ing knowledge of our globe and its gravity anomalies, the descriptive 
geographical parameters are subject to change, particularly the number 
and values of the oblateness coefficients, but also -- to a lesser 
degree -- the gravitational parameter GM©, and the geometrical param- 
eters flattening f, and mean equatorial radius Rg. Urgently needed 
standardization of the principal parameters was accomplished by NASA 
in 1963, but new and additional data have been determined continuously 
ever since and are being determined from an ever-increasing population 
of artificial satellites in precis ion- tracked (Baker-Nunn) orbits of 
all inclinations. 

To relate these new observations to the geoid without having to 
accept the inconsistencies of the second- and third-order Clairaut 
extensions, a new theory was developed recently at MSFC by H. Krause^®’ 
which further increased the bulk of existing geographic parameters and 
models . 

It is the objective of the following analysis to define and to com- 
pare seven different geographic models, one of which is the current 
(but outdated) NASA Standard. These models are based on four different 
bets of oblateness coefficients, all of which are applicable to present 
practical work, and also on the fact that for each set of oblateness con- 
ditions two body geometries can be defined; i.e., a spheroid composed of 
superimposed spherical harmonics and a regular ellipsoid approximating 
this spheroid. 

In the first part of the analysis, the general background and the 
nature of the problem is considered in some detail, followed by the 
derivation and description of the geophysical theories. In the third 
part, the models are compared with each other. A number of salient 
conclusions, forming the fourth part, can then be drawn. 
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II. GENERAL BACKGROUND AND DISCUSSION OF THE PROBLEM 


A. The Potential Function of the Earth 


In the study of geocentric motion such as that encountered on rocket 
trajectories and satellite orbits, among the most important constants 
necessary for describing the environment are those pertaining to the 
gravity field acting on the moving body, and to the shape or figure of 
the gravitating body, the earth. 

Theories attempting to describe the earth's gravity field are, by 
necessity, fundamentally based on Newton's Universal Law of Gravitation, 
which is rigorously accurate for a central force field as produced by a 
perfectly homogeneous and perfectly spherical body. The gravitational 
potential of such a Newtonian force field is 


U = - . (1) 

The potential function of the earth U 0 , or the geopotential, is 
defined as the integral of the gravity forces over the entire field, 
which are the resultant of the gravitational forces due to Newtonian 
attraction and the centrifugal force due to the rotation of the earth 
and its atmosphere. Thus, the geopotential is the sum of the potential 
of the gravitational field and the potential of the centrifugal force. 


U 0 = ■ r (1 + ^) + “I cos^ 0 . 


( 2 ) 


If the rotational term is omitted, equation (2) expresses the 
gravitational potential which is appropriate to inertial coordinates, 
as useful for exo-atmospheric (free-flight) flight phases. For geodesy, 
as well as for all other applications requiring a rotating earth coordin 
ate system, the, term is retained. 
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In equation (2), the modifier i]/ accounts for the fact that, contrary 
to the basic assumptions of the Newtonian potential, the earth's internal 
mass distribution is not entirely homogeneous and its shape is not 
spherical, but oblate. Also, its surface is covered with irregularities, 
such as continental highlands, depressions, mountain ranges, ocean deeps, 
valleys, etc., which give rise to gravitational anomalies and require an 
additive perturbative term to complete the geopotential. 

Since the potential function of the earth is a solution of the 
Laplace equation (V®U^ * 0), it is a harmonic function and can be repre- 
sented analytically as the sum of a series of two-dimensional (spherical) 
harmonics in a manner similar to the one-dimensional Fourier series as 
representation of a nonanalytical function on a circle. The coefficients 
of the terms of this infinite (converging) series are polynomials of the 
general Legendre type, while the terms themselves are the so-called sur- 
face harmonics, the introduction of which is due to Legendre and Laplace. 
If the spherical harmonic is a function of two variables, such as lati- 
tude and longitude (rather than only one, the latitude), it involves the 
so-called Associated Legendre Polynomials. 

The acceleration of gravity g of the earth is defined as V U and 
can therefore also be expressed in spherical harmonics. 

The perturbative parameter of the gravitational potential in 
equation (2) can be expanded, as described, in an infinite series of 
spherical harmonics. 


00 
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n=l 


i (-^T >:»■ » 

m=0 


C cos mA + S sin mA , 
n,m n,m * 


( 3 ) 


where r is the distance from the center of the earth, is the earth's 
mean equatorial radius, 0 is the latitude, A is the longitude, and 

Sn m numerical coefficients, and the are the associated Legendre 
polynomials, Introducting equation (3) in equation (2), one obtains 
the general formula for the earth potential, as recommended by the lAU 
( 1962 )^’^. 
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( 4 ) 



P”(sin 


cos roA + sin mA) 
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cj r 


cos^0. 


B. The Spherical Harmonics 

It is seen that the second (aspherical) part of equation (4) con- 
sists of an infinite series of harmonics which -- in a rigorous sense — 
describes the perturbative potential accurately only if an infinite 
number of terns are considered. Taking the terms by families according 
to the indices, it is furthermore seen that for the first existing 
family, where m = 0, 1 < n < », the equation reduces to 
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i-r) ■ ® • <=n.' 


n=l 

for the purely gravitational potential. 


(5) 


In this special case, P^Csin 0) are the (simple) Legendre poly- 
nomials of argument sin 0, usually referred to as "zonal” harmonics . 
Obviously, the potential Um»o is independent of the longitude A; hence, 
it is constant along (latitude) parallels. This means that the undula- 
tions of the hairmonics on the sphere are zero along parallels and become 
alternately positive and negative between these n parallels in a dis- 
tance jt. Thus, equation (5) refers to an axially symmetric earth which 
is fluted along the parallels and where the symmetry axis is the rota- 
tional axis through the poles. 

As an alternate for Cn,0» it ^s been recommended by the lAU^ to 

use 



( 6 ) 
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Thus, for the body-of-revolution case, the purely gravitational 
potential becomes 


U 


m =»0 




P^(sin 0) . 


( 7 ) 


The true figure of the earth does not exactly exhibit rotational 
symmetry, as assumed by setting m ■ 0 in equation (4). The principal 
axes of inertia of the true earth do not exactly coincide with the 
rotational axis system. Thus, the angular momentum of the earth is 
not pointing exactly in the same direction as the angular velocity. 
Consequently, products of inertia do exist, which are represented by 
higher harmonics. By letting m assume integral values from 1 to n, 
terms are added to the second part of equation (5), which include the 
longitude A as a second variable. Consequently, the associated 
Legendre polynomials P^(sin 0) are also required, as well as coef- 
ficients Jn,m* These terms are referred to as tesseral harmonics , with 
the exception of the two terms in equation (4) resulting from letting 
m = n, so that 


C p"^ cos nA and S P*^ s in nA 
n,n n n,n n 

appear. These two terms are called " sectorial" harmonics . 

The tesseral harmonics (0 < m < n) have zeros both along meridians 
and along parallels; thus, the undulations represented by them become 
zero on the reference spheroid simultaneously on a number of meridians 
and parallels. They undulate in the form of a network across the spheroid 
similar to a chessboard composed of fields that are alternatingly positive 
and negative. Thus, for example, the spherical harmonic represented by 
the function Ps (with the coefficient Js, 4 *) forms a network of five 


* The notation used commonly in the literature for the tesseral and 
sectorial J-coeff icients is Jnm- Thus, the first sectorial harmonics coef 
ficient would be Jgsj first tesseral coefficient J 2 . 2 , etc. While this 

practice represented an ambiguity in notation, with respect to the zonal 
harmonics, no great confusion could be caused as long as the zonal har- 
monics used in practical work did not exceed the tenth order. With pre- 
sent analyses of the higher zonal harmonics already going to n ■ 21 , 
m = 0 ,^ however, notation of the mentioned type would easily confuse 
tesseral with zonal harmonics. Thus, the zonal J 14 could easily be inter- 
preted as the tesseral n = 1, m = 4. Since powers of are sometimes 
encountered, the old practice of writing is also not desirable. Con- 
sequently, it is suggested that tesseral and sectorial coefficients are 
indexed with a dividing comma, i.e., Jn,m» while the zonal notation 
be retained. The notation P^ for the associated Legendre polynomials 
should also remain unchanged. 


meridians (four zeros in 180°) and (6 - 4) = 2 parallels (not equator) 
on the sphere, cutting it into 30 (alternatingly raised and depressed) 
pieces. 

In the special case of the two families of sectorial harmonics, the 
undulations of the function become zero on the reference spheroid for 
cos 0 = ± 1 (on the poles) and for cos nA = 0, i.e., along n meridians 
which are symmetrical to n planes through the rotational axis, thus 
cutting the spheroid in sectors which are alternately concave and convex. 
J4 4, for example, would represent undulations featuring four meridians 
and eight sectors. 


C. The Harmonics Coefficients 


The coefficients of the tesseral harmonics of the earth, especially 
the higher ones, are not really known. The sectorial coefficients J2, 2 
and J4j4 have been determined tentatively during the past 6 years and 
are sometimes used in high precision orbit determination programs®. For 
J2, 2, which has two zeros along a parallel (four sectors) and thus deter- 
mines the ellipticity of the equator and the parallels, Kaula® gives 


J2, 2 “(1.80 ± 0.1) X 10"®. 


Of the tesseral harmonic perturbations of satellite orbits, only 
the sectorial J2 2 1-s large enough to be of practical concern in most 
orbit analyses; specifically, for example, this harmonic affects the 
supplemental energy requirements of 24-hour satellites. The sectorial 
harmonic due to J2, 2 depicted in figure 1. 

Because of the uncertainty in and relative insignificance of the 
tesseral harmonics coefficients, they are usually neglected in the 
potential function which therefore assumes the commonly used form 
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n=l 



P„(sin 0 ) 


( 8 ) 


which was first suggested by Brouwer^ and adopted by the lAU in 1961. 

If so desired, tesseral and sectorial terms with the appropriate assoc- 
iated Legendre polynomials can be added to the second part of the equa- 
tion, as required by the inclusion of the longitude A in the analysis, 
as well as the rotational potential (l/2)r^oj^cos^0 to the whole equa- 
tion to transform the gravitational potential in inertial coordinates 
for a free body to the gravity potential for a body attached to the 
earth. 
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( meters ) 



DUE TO 499 ) at equator 






of the (theoretically infinitely many) zonal harmonics coefficients 
Jjj, only a few are known with adequate accuracy. If an equatorial geo- 
centric coordinate system is assumed as a reference system such that its 
origin coincides with the center of mass, and if the latitude 0 is meas- 
ured relative to this center, the first term in the sum of equation (8), 
which is proportional to 1/r, vanishes since the coefficient J^, 
multiplied with Rg, would represent the distance between coordinate 
origin and center of mass; thus, = 0. 

Of the other zonal functions, the second order harmonic is by far 
the largest, since it refers to the mass distribution due to the earth 
oblateness, thus defining the departure of the shape of an earth ellip- 
soid from that of an ideal sphere, while J4 being three orders of magni- 
tude smaller than J2, and all following even J's account for deviations 
in hemispherical mass distribution symmetrical about the equator plane, 
i.e., for deviations of the geoid from the ellipsoid. Each even J of 
increasing order, while maintaining symmetry about the equator plane, 
"narrows down" the figure into an increasingly complex shape, while the 
sum of subsequent even j's describe the deviation of the geoid from the 
preceding shape. The even J*s, however, all represent spheroid shapes 
with identical northern and southern "hemispheres." 

The asymmetry of the northern and southern "hemispheres" of the 
geoid, also referred to as "pear" shape, is taken into account by the 
odd surface harmonics, which contain odd powers of sin 0 , especially by 
the principal coefficient, J3. Higher-odd-order j's are at least one 
order of magnitude smaller than this coefficient. 

The first zonal harmonic is shown in figure 2 "^. Figures 3 and 4 
depict the zonal harmonics due to J3 and J4, and J5 and Jg, respectively. 
Of particular interest is the ovoid form of J3, which is responsible for 
the "pear" shape of the earth. 


D. The Theorem of Clairaut 


In the past, higher order zonal harmonics, and especially the odd 
harmonics, have been neglected in the geopotential. This was mainly 
due to the methods of determining the coefficients used before the advent 
of the artificial earth satellites. Also, since the amplitudes of high 
harmonics -- according to the quantity l/r"^ in their terms -- decay very 
rapidly with altitude above the surface, and their effect on bodies in 
the earth's external gravitation field dwindles rapidly with increasing 
order and/or altitude, coefficients above fourth order, as well as the 
odd term J3, were usually neglected. Before the advent of the artificial 
satellites, the moon was the only moving body in the earth's potential, 
but at its distance only the second harmonic is still discernible in its 
motion. 


* The figure is actually the second term of equation ( 34 ), 
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FIG. 2. SECOND SPHERICAL HARMONICS TERM 
[-14266 Pt (sin ^)] IN EARTH 6E0I0 
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FIG. 4 . FIFTH AND SIXTH ORDER SPHERICAL HARMONICS 
TERMS [1.3 P5 (sin - 4 Pg(sin^)] IN EARTH 6EOID 



Determinations of the J-coeff icients Jg and J 4 were fundamentally 
based on geodetic surveys, gravity measurements, observations of the 
moon, and Clairaut's theorem. The latter deals with the form of a sur- 
face which encloses all the matter of a rotating body with various dens- 
ity distributions, which constitutes an equipotential surface (surface 
of constant potential). Clairaut's formula (1743) is 

J = f - I m, (9) 


where J is the second harmonics coefficient, m is the centrifugal factor 
at the equator co^R/gg (with = centrifugal acceleration at equator, 
gg = gravitational acceleration at equator), and f is the flattening of 
the ellipsoid of revolution, as assumed here. 

Obviously, Clairaut's formula relates the amplitude of the second 
harmonic in gravity to that of the corresponding harmonic in the radius 
vector of the equipotential surface. 

Clairaut's formula, because of the approximations made in deriving 
it, is not accurate enough for practical use. A second order extension 
of the theorem has been derived by Darwin (1899), DeSitter (1924), 
Jeffreys (1954), Sullen (1946), Cook (1959), et al., resulting in equa- 
tions involving higher orders of the flattening, such as the third-order 
equation given by Kaula^, 
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f2 + 0(f3) 


( 10 ) 


With the use of Clairaut's theorem, as well as its extensions by 
the rtientioned authors, the earth flattening -- before the artificial 
satellites -- could be derived separately in two ways; from a harmonic 
analysis of observed gravity values (where the measured gravity anomalies 
are used to compute the absolute deflections of the vertical and, thus, 
the true undulations of the geopotential surface; also, to derive the 
detailed gravity field) and from astronomical observations of the con- 
stant of precession and the theory of the moon's motion (DeSitter, 1915), 
leading to Jg. Higher j' s could then be computed from the higher-order 
extensions of Clairaut's theorem, using the derived flattening, but when 
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the first artificial satellites were flown, it was found that higher* 
order harmonic coefficients determined with this method did not agree 
with the J*s actually observed from the orbital perturbations of the 
satellites. The reason for this serious discrepancy is simply the fact 
that for the higher-order terms, due to the irregular density distribu- 
tion within the earth, the basic assumption of an ellipsoid of revolu- 
tion as the figure of the earth in the Clairaut formula and all 
higher-order extensions is an increasingly bad approximation of the 
geoid . 


E. The Geoid 


Distinction is made between the geodetic geoid, which depends on 
the surface irregularities of the earth, and the geophysical geoid, which 
is governed by the internal irregular mass distribution of the earth. 

The geodetic geoid is the true figure of the earth. While the actual 
earth surface with its irregularities is not an equipotential surface, the 
geophysical geoid, however, is required to be an equipotential level, 
identical with mean sea level and its continuation under the continents. 
The geoid, thus, is an extremely complex shape, requiring a very large 
number of higher-order harmonics for its description. 

The reason for the use of the geoid in geodesy is the following: 

In principle, the geometrical form of the earth’s surface can be found 
independent of observations of the potential and gravity, by geodetic 
triangulation methods, and then, given the values of these quantities 
on the surface, Laplace’s equation can be solved for all the exterior 
space. Because of the extreme complexity of the true boundary, this 
would be an almost hopeless task. It is therefore universal practice 
to refer all geodetic observations to the gravitational equipotential 
surfaces and to determine the form ,of ^ these from geodetic observations. 
Because the sea-level equipotential surface is an internal surface on 
land, the gravity field is not computed for the actual earth but for a 
model earth which is related to the actual earth and which is bounded 
by an equipotential surface. The field so computed will not agree with 
that of the actual earth throughout space, but the model may be chosen 
so that it agrees where observations can be made. This is the geoid. 

The geodetic geoid must be computed point by point and cannot be 
given* by' few parameters, since it is not a mathematical surface but 
depends on the irregular distribution of visible and invisible masses 
near the earth’s surface. For geodetic surveys, however, this point-by- 
point calculation would be impractical. Therefore, reference spheroids 
of revolution have been used instead. The reference spheroid is assumed 
to be an equipotential surface of the same volume and flattening as the 
geoid , 
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The history of geodesy has seen the following principal reference 
spheroids (all ellipsoids): 

Germany (and several European states) 1841 Bessel 

England 1880 Clarke 

U.S.S.R 1938. . . , .Krassowski 

U.S 1866 Clarke 

U.S 1910 Hayford 

U.S 1963 Kaula/Fischer 


F. Determination of the Harmonics Coefficients 


As has been mentioned, Jg in former times was determined from the 
moon's motion. For the lower amplitudes of the higher order harmonics, 
however, these effects are not discernible because r becomes too large. 
However, when the first artificial satellites were orbited, it was found 
that the gravitational anomalies expressed by the higher-order harmonics 
characterize the potential field at their altitudes. Thus, the existence 
and magnitude of J 3 (i.e., the "pear" shape of the earth) was first dis- 
covered from satellite 1958^2 (Vanguard-1), launched in March 1958®. 

The asphericity of the earth, represented by the spherical harmonics, 
causes constant Secular and periodic perturbations on a satellite orbit. 

By observing the secular and periodic perturbations, the zonal harmonics 
coefficients can b© determined. Since the secular changes, especially 
the orbital precession (i.e., the regressive rotation of the nodal line) 
and the advance of the perigee (i.e., the rotation of the apsidal line), 
depend primarily on the even-numbered Jjj, while the long-period oscil- 
lations in four of the five orbital elements (eccentricity, inclination, 
longitude of ascending node (right ascension), and argument of perigee) 
are caused primarily by the odd harmonics®, the determination of the 
even is usually done from secular changes, while the odd Jn are com- 
puted from periodic perturbations. An orbital perturbation suffered by 
one particular satellite yields one linear equation between the harmonics. 
For example, if p is the amplitude of the observed oscillation in eccen- 
tricity, we obtain an equation for the coefficients of the odd harmonics, 
^ 3 > Jsj J 7 . . . of the form‘d 


A 3 J 3 + A 5 J 5 I" A 7 VJ 7 ... ~ 3, 


( 11 ) 
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where the A's are constants for a particular orbit and depend mainly on 
the orbital inclination i and semi-major axis a. By using several dif- 
ferent satellites, several equations of the form (11) can be obtained. 
Therefore, if a large number of satellites (widely distributed over as 
many inclinations as possible) is available for analysis, an equally 
large number of can be determined. 

In solving the simultaneous equations (11), the main error of the 
analysis is introduced by the number of to be determined remaining 
finite and all terms higher than the highest desired J being assumed 
negligible. In fact, however, these higher coefficients are not neces- 
sarily negligible, and since their effect is - by necessity - lumped 
together in the computed Jn, in any one analysis of a given set of 
satellites there are always several different solutions to the J's, 
depending on the number of j's considered in the algebraic equations. 
Thus, for example, King-Hele, et al.'^, recently computed odd zonal har- 
monics up to J21 from eccentricity oscillations of 14 satellites, giving 
two different sets, i.e., one with J3 to J]_5, and the other with J3 to 

Jai. 


Also, the equations are sometimes solved by several different 
methods to obtain more confidence in the results. For example, the 
equations may be solved for three j's at a time, assuming the higher 
j's to be zero, in fours and fives, and also with least-squares and 
minimum-residual methods. The coefficients given by King-Hele'^ in the 
larger set are 




J 3 

= 

(-2.50 

+ 

0.01) 

X 

10"® 



J 5 

= 

(-0.26 

± 

0.01) 

X 

10-® 



J7 

= 

(-0.40 

± 

0.02) 

X 

10"® 



J 9 

= 

( 0 

± 

0.06) 

X 

10-® 



J 11 

= 

(-0.27 

± 

0.06) 

X 

10"® 

>. 

(12) 

J 13 

= 

(+0.36 

± 

0.08) 

X 

10-® 



J 15 


(-0.65 

± 

0.10) 

X 

10"® 



J 17 

= 

(+0.30 


0.08) 

X 

10-®. 



J 19 

= 

( 0 

± 

0.11) 

X 

10-® 



J 21 

= 

(+0.58 

± 

0.11) 

X 

10"® 

J 



The latest most 

: complete 

set of 

zonal harmonics, 

both even and odd 


coefficients, has been evaluated by Kozai. In 1963, Kozai^° derived a 
set of values for eight coefficients of zonal harmonics (up to and 


16 



and Including Jg) from the (then) available observations of secular 
motions of the node and the perigee and the amplitudes of long-periodic 
terms of artificial satellites. Later, in 1964, from precisely reduced 
Baker-Nunn observations of nine high-inclination satellites (28® - 95®), 
he produced*^ a new set of thirteen coefficients (up to Jx 4 )j 


J 2 

S 

(1082.645 

J 4 

= 

(-1.649 

Je 

= 

(+0.646 

Ja 

= 

(-0.270 

J 10 

= 

(-0.054 

J 12 

= 

(-0.357 

Jl4 

= 

(+0.179 


± 0.006) X 10"® 
± 0.016) X 10"® 

± 0.030)x 10"® 
± 0.050)x 10"® 
± 0.050)x 10"® 
± 0.047) X 10*® 
± 0.063) X 10"® 


J3 =(-2.546 
J5 =(-0.210 
J-7 =(-0.333 

Jg =(-0.053 
Jii =(+0.302 
Ji3 =(-0.114 


± 0.020) X 10"® 

± 0.025)X 10"® 

± 0.039) X 10"® 

± 0.060) X 10-® > 
± 0.035) x 10"® 

± 0.084) X 10"® 

J 


(13) 


The values of Kozai's even coefficients have been compared recently^^ 
with two other sets of King-Hele/Cook and Smith. The agreement is excel- 
lent for all inclinations greater than 28®, which was the smallest incli- 
nation of available satellite orbits. It is therefore concluded that, 
as long as there are no satellites at inclinations between 10® and 25®, 
the evaluation of the even harmonics in the potential carried beyond the 
present status does not seem worthwhile. For the odd harmonics, this is 
not the case, since their primary effects, the long-periodic changes 
(e.g., the amplitude of the oscillation in eccentricity) decrease to 
zero as' the inclination tends to zero. • For this reason, King-Hele felt 
justified to publish the odd zonal harmonics to Jsi* as given above, with 
an indication that J 23 and J 25 are small. 


G. Krause's Theory of the Geoid Surface 
If Kozai's value for J 2 , 


J2 =(1082.645 ± 0.006) x 10‘®, 


is introduced in a relation of a third-order Clairaut theory^. 



1 


- 2 t + + ^ + 0(£3)j. 

(14) 
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a value of f = 1/298.254 results for the flattening of the reference 
earth ellipsoid. Using this f to compute the theoretical value of J 4 
for the reference ellipsoid assumed to be in hydrostatic equilibrium, 
one obtains -2.350 x 10 “ ®, as compared to the observed value of 
-1.649 X 10'®. This discrepancy illustrates the increasing inaccuracy 
of Clairaut's (ellipsoid) theorem and its extension, if applied to 
higher-order harmonics. 

To find a more accurate relation between the radius vector of the 
earth and its gravity field, thereby cleaning up the inconsistency in 
the older theories, Krause recently developed a nexvr theory^®* which 
essentially introduces the radius vector of the geoid . rather than of 
an ellipsoid, consistently expressed to the same accuracy as the poten- 
tial field it is related to. Since the new theory is kept general, it 
can handle all higher order zonal harmonics up to Jn with the same 
internal consistency, which the older theories, based on Clairaut, Darwin, 
Helmert, et al., are unable to do. Using Krause's theory with Kozai's 
thirteen coefficients and introducing the radius vector of the geoid 
expanded in spherical hamonics up to 3^4 in the similarly expanded geo- 
potential, a new flattening of the earth can be computed, f = 1/298.1840*, 
which differs both from the old second-order value given above and from 
the flattening of the Kaula/Fischer ellipsoid used currently as a world 
datum, for example, in trajectory computations (1/298.30). It is 
believed that the value derived with Krause's theory is more consistent 
with the observed spherical harmonics coefficients. 

In his theory^®’ Krause has made use of the fact that the 
generalized (equipotential) surface of the geoid, similarly to the geo- 
potential, can be represented alternatively by two series expansions, 
one consisting of spherical harmonics, the other of powers of the sine 
of the latitude. When introduced in the geopotential, the chosen series 
must be of the same degree as the potential function used. The author 
shows that the coefficients of the terms of the series (Ajj in the 
spherical harmonics series, in the power series) are related to the 
known oblateness coefficients Jn of the geopotential equation U and can 
be computed from these. Insertion of the radius equation thus obtained 
in the equation of the total gravity acceleration in the direction of 
the normal yields the gravitational acceleration of the geoid as an 
expansion of surface harmonics, based on a consistent geoid geometry. 
Because of the relation between the Ajj (or Bn) coefficients and the 
oblateness coefficients Jn, and assuming the surface to be of constant 
potential, the oblateness of the northern and southern "hemispheres" and. 


"The value given in references 13 and 14 has been slightly corrected by 
the present writer. 
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thus, the mean meridional flattening f can be computed from the coef- 
ficients, relating now to the true geoid and the full set of observed J„, 
rather than to a reference ellipsoid, as in Clalraut's theory. 

For the geoid radius, expanded in spherical harmonics. 


00 



- AjjPjjCsin 0 ) 

n=l 


+ 3A2,2 CO8^0 cos 


2(A - V. 


(15) 


in which the first sectorial harmonic, with the associated Legendre poly- 
nomial P 2 (sin 0) exactly given by 3 cos^0, is taken into account to express 
the ellipticity of the earth equator, Krause finds for the A-coeffi- 
cients^^* 


Aq ~ 1 


1 . 

X 



l>3«5...(2v - 1) 
2»4*6. . . 2v 


2V 








2 “ 


Wg ■ Bs) 
(16) 



1 

J2 + 3 Wg 




- Bs) 


(17) 



J4 









(18) 


A 


K 


1 

X 


J 


K 


(K 0,2,4) 


(19) 


where 





= centrifugal factor at equator, 


( 20 ) 


Bp - 


v=0 


(- 1 ) 


V 1*3»5. . . (2v + 3) 


2^ v! 


1 ~ 

>7 "o 

2V^2 I e 


( 21 ) 
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and 


X = 1 


00 



n=l 


3-5-7,. . (2n + n 
2-4.6... 2n ■'^sn 





( 22 ) 


If the first variation of U with longitude, as represented espec- 
ially by the first sectorial term, is taken into account as it was in 
equation (15), where 3 cos^0 is used as an exact substitute of Pg(sin 0) , 
the geopotential of the earth becomes in inertial coordinates 


U 



/ R sn / R X ^ 

—^J Pjj(sin 0) + \—^J cos20 cos 2(A - A^) 


+ 


2 


U) 

e 



(23) 


Since the three components (in polar coordinates) of the gravity 
are defined as 


= ^ 

®r cir 

= 1 M 

®A r cos 0 cl A 

- i ^ 

^0 r b0 


(24) 


while the total gravity acceleration in the direction of the normal is 


8 




(25) 
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the gravity. Including the sectorial term, can be given by determining 
the potential derivatives and expanding the sum of their squares according 
to equation (25) In a series approximating the square-root of equation 
(25). The result ls^®>^^ 


g 



n-2 



• P^(sln 0) 



co8^0 cos 2(A - 



e 



CO8^0 


(26) 


where only the derivative of P 2 (sln 0) with respect to 0 has been 
Included (expressed as 

dP (sin 0) 

P 2 (sln 0) » * 3 sin 0 cos 0), 


while the higher order derivatives have been neglected. 


III. ANALYSIS OF MODELS 


A. Model 1 (Ja-Only) 


For computation of trajectories close to earth, featuring small 
central range angles, the effects of the equatorial bulge and other 
Irregularities can usually be considered negligible, primarily because 
the flight times involved are relatively short. In this model, it is 
therefore assumed that the earth is an ellipsoid of revolution, sym- 
metrical around all three axes. Except for the second harmonic, all 
higher order effects of asphericity (which are three or four orders of 
magnitude smaller than are taken as zero. 
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With these assumptions, the gravity equation (26) reduces to 



3J2 



P2<sin p) - 



cos^0 



sin^0 008^0 


( 27 ) 


It would be more accurate to compute the gravity directly from 
the potential of the ellipsoid, rather than from the series approximation, 
equation (27). The task of determining the potential of a homogeneous 
ellipsoid was a very formidable one and remained unsolved till the end 
of the 18th century, despite attempts at a solution by Newton, Maclaurin, 
D'Alembert and Lagrange. Finally, Laplace succeeded in finding the 
general solution, followed later by Gauss and others. The general form 
of the potential of a homogeneous, triaxial ellipsoid for an exterior 
point involves an elliptic integral of the first kind; hence, the com- 
ponents of the gravity force of the ellipsoid also contain elliptic 
Integrals. For the present case of the ellipsoid of revolution, the 
equations become considerably simpler; the elliptic integrals reduce to 
logarithmic-cyclometric integrals. The investigations of the above men- 
tioned authors generally refer to homogeneous bodies. Clairaut was the 
first who investigated inhomogeneous bodies. 

A third-order relation of Clairaut' s theory is given by Kaula^; 


T _ 2 1 1 

02 2 ^(^ “ 2 ^ m 


-I 


m 




+ ^ m^ 
4 


ii f2 

49 


+ 0(f3) 




(28) 


Kaula gives also 

N 

GM^= 3.986032 x 10^° cm^ sec"^ 

= 398603.2 km^ sec"^ 

w “ 0 . 729 , 211,585 x lO"^ sec'^ ) . 

Rg = 6378165.0 m 

g = 978.0300 cm.sec*^ 
e J 


( 29 ) 
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Using these values, the centrifugal factor at the equator, m, 
becomes 


• R 

m « ^ = 0.003,467,773,255, (30) 

®e 

which is in agreement with reference 6. 

For Model 1, the value for Jg used is the first one of the set given 
recently by Kozai^^, equation (13): 


Jg = 1082.645 X 10- ®. 


(31) 


In a strict sense, this Jg-value is not completely consistent with 
the above given value of gg from Kaula^. The equatorial gravity of the 
oblate spheroid may be computed from 


g 


e 


0 • 



9 


where the factor a accounts for two effects. First, it includes a cor- 
rection to M0 because the mass of the earth's atmosphere, included in 
M0,does not contribute to the gravitational acceleration at the sur- 
face, and second, it modifies the succeeding term, valid for a sphere 

only, for the oblateness effects. In the Kaula value of equation (29), 

0 by necessity refers to the older Jg and probably also to the older J4. 
However, the inconsistency is not felt to be of practical consequence 
for the present purpose. For the Jg of equation (31), 0 would become 

0.9981616959 (by solving equation (22) and computing 0 = X(1 - A), where 

A = Matm/M0* 0.8594 x 10"^®)^^, so that the above equation results in 
gg = 978.026 gal, 4 mgal lower than Kaula* s value, which is presently 
accepted as standard^®. 

It must be kept in mind that the assumption of all higher being 
zero introduces, in the strict sense, an error in the model which is not 
taken into account by the above Jg-coefficient, since its particular 
value has been determined from simultaneous algebraic equations involving 
12 additional coefficients of higher order, which were different from 
zero. However, it is assumed that the uncertainty caused by this fact 
will be of no concern for the present purposes. This assumption is 
justified by the results of this analysis. 
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Introducing the Jg-value in equation (28), one obtains the oblate- 


ness 


f = 0.0033528465 = 1 : 298.254, 


which is consistent with the assumed geopotential. 

With the value for Rg given in equation (29), the polar earth 
radius becomes 


Rp = RgCi - f) (32) 

R = 6356780.0 m 
P 

R = 6378165.0 m 
e 


and the general radius at latitude 0 


R R 

R = - ' ■ ■ ^ P (33) 

n/r cos^0~^^ s in ^0 

p e 


Since the assumption of an ellipsoid is not fully consistent with 
the gravitation assumed (based on Jg), a second body has been deter- 
mined by using the first spherical harmonic in a truncated expression 
of the type (15). By solving equations (16), (17), (20), (21) and (22), 
the coefficients are determined: 


Aq = 0.9988810301 

Ai = 0 

A 2 = 0.2236813952 x lO"^ 

A3 = 0 

A 4 = -0.1501014034 X 10-s. 
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Neglecting the coefficient A 4 , the radius of the spheroid id then 
determined from 


Ri = 6371028.026 - 14266.76846 • P 2 (sin 0 ) , (34) 


where PaCsin j2) is the second -degree Legendre polynomial. 

The resulting spheroid is depicted in figure 5. 

The total gravity (equation (27)) and the surface shape of the 
ellipsoidal Model 1 (equation (33)), as well as the surface of the body 
formed by the Ja-harmonic (equation (34)) are given in table 1. 


B. Model 2 (NASA Standard) 

For trajectories of longer duration, such as those of upper stages, 
as well as for near-earth orbits, the geopotential exhibits anomalies 
such that its effect deviates from that of an attracting ellipsoid of 
revolution. For almost all practical purposes, however, it suffices to 
add only the two next higher-order zonal harmonics to the potential and 
to consider the gravity anomalies averaged over all longitudes, thus 
retaining the rotational symmetry around the polar axis and rendering 
tesseral and sectorial terms unnecessary. 

The second model is represented by the geopotential and the astro- 
geodetic world datum as adopted by the Ad Hoc NASA Standard Constants 
Committee^® in 1963 for Project Apollo and other NASA programs. The 
standard is based on the Kaula/Fischer Ellipsoid of 1963. Irene Fischer 
of the U. S. Army Map Service, in 1960, published a world ellipsoid 
based on an imposed flattening f =■ 1/298.30 (which at that time appeared 
to be the best available value from early artificial satellites), by 
determining the ellipsoid of revolution which best fit the geoid con- 
tours derived from astrogeodetic measurements^®. Kaula later^, in 1963, 
published a new value for the equatorial radius. 


R * 6378165.0 ± 25.0 m, 
e 


(35) 


which was a compromise between Fischer's value^® and a 1961 value by 
Kaula^^. The NASA Committee meeting at Goddard Space Flight Center in 
May 1963 adopted this value, along with the following constants^®: 
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FIG. 5. HEIGHT OF J 2 - SPHEROID (MODEL 1) OVER REFERENCE 
ELLIPSOID OF FLATTENING 1/298.222 (MODEL 3) 


26 








TABLE 1 


Model 1 


Latitude 
North or South 
(deg) 

Radius 

(Ellipsoid) 

(m) 

Radius 

(Spheroid) 

(m) 

Height of Sphe- 
roid over 
Ellipsoid (m) 

Gravity 

(gal) 

90 (poles) 

6356780.0 

6356761.3 

-18.7 

983.2049 

85 

6356941.6 

6356923.8 

-17.8 

983.1663 

80 

6357421.7 

6357406.6 

-15.1 

983.0516 

75 

6358205.8 

6358194.8 

-11.0 

982.8641 

70 

6359270.5 

6359264.6 

- 5.9 

982.6093 

65 

6360583.7 

6360583.5 

- 0.2 

982.2944 

60 

6362106.1 

6362111.3 

5.2 

981.9286 

55 

6363791.7 

6363801.7 

10.0 

981.5228 

50 

6365589.7 

6365603.3 

13.6 

981.0890 

45 

6367445.6 

6367461.3 

15.7 

980.6401 

40 

6369303.1 

6369319.4 

16.3 

980.1897 

35 

6371105.7 

6371121.0 

15.3 

979.7517 

30 

6372798.5 

6372811.4 

12.9 

979.3394 

25 

6374329.7 

6374339.2 

9.5 

978.9656 

20 

6375652.3 

6375658.1 

5.8 

978.6422 

15 

6376725.7 

6376727.9 

2.2 

978.3794 

10 

6377517.0 

6377516.1 

- 0.9 

978.1854 

5 

6378001.7 

6377998.9 

- 2.8 

978.0664 

0 (equator) 

6378165.0 

6378161.4 

- 3.6 

978.0300 
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GM© = 398603.2 km sec"^ 
f * 1 : 298.30 
gg = 978.030 gal 

J = 1.62345 (±0.00030) x 10-3 

. 

H = -0.575 (±0.025) X 10’® / 

D = 0.7875 (±0.0875) x 10"^ 

R « 6356783.3 m 

P 

u =» 3461.4143 X 10 ”® sec”^ 
e 


( 36 ) 


where the J, H, D coefficients apply to an alternate (older) form of the 
geopotential, and the value 


u 

e 


w^R® 

e 

GM 


is fully consistent with the assumed parameters. 

In the modern notation, the oblateness coefficients become 


Ja = I J = 1082.30 X 10”® 
J 3 = I H = -2.30 X 10”® 

J 4 =--^D = -1.80 X 10”® 


(37) 


It is pointed out that the above set of constants is consistent 
within itself for most practical purposes. However, since the geometrical 
figure assumed is still an ellipsoid of revolution (equal northern and 
southern "hemispheres"), derived from J 2 through a third-order relation 
of Clairaut's theory (equation (28))^, it is at best only an approximation 
and contains theoretically an inconsistency. 
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Table 2 lists the gravity of Model 2 (second column), computed 
from equation (26) without the longitude-dependent term, with J 2 , J 3 , J 4 
as given by equation (37), and in the third column the local radius of 
the associated Kaula/Fischer reference ellipsoid, computed from equation 
(33) with the data given in (35) and (36). 

As an illustration of the inherent inconsistency of the NASA Standard 
model, it is seen that the gravity field at the North Pole is higher 
by 18.4 mgal than at the South Pole, while the radial distance computed 
from the model is unable to account for this "pear" shape. 

To compare the deviation of the ellipsoid from the (here assumed) 
geoid, the shape of the geoid must be computed. The new theory by 
Krause^^’^^, discussed in Section II, provides a relatively simple method 
for doing this. To write the geoid radius in an expansion of the second, 
third and fourth surface harmonics, the coefficients An in equation (15) 
must be determined with J 2 , J 3 , J 4 (equation (37)), from equations (16) 
through ( 22 ). 


One obtains 

Aq = 0.9988805375 '' 

Aj^ 0 

A 2 = 0.2236440156 x lO'^ > , 
A 3 = -0.2304227311 x 10"^ 

A 4 = -0.3313107217 X 10'^ ^ 


(38) 


with which the geoid surface of Model 2 can be expressed by 


Ra = 6371024.88 - 14264.38 • P 2 (sin 0 ) + 14.69 • P 3 (sin 0 ) 

+ 21.13 • P 4 (sin 0 ) (m) (39) 

where the Pjj(sin 0 ) are the Legendre polynomials of argument sin 0, 

Column 4 of table 2 lists the results of equation (39), while the 
deviations of the ellipsoid from the geoid are given in column 5. It is 
seen that the NASA Standard Model results in an earth radius which is 13 m 
below the geoid at the North Pole, 16.4 m above the geoid at the South 
Pole, and between zero and 12 m above or below it in the intermediate 
latitudes . 
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TABLE 2 
Model 2 


Latitude 

(degrees) 

Gravity 

(gal) 

’ 

Rad lus^ 
Ellipsoid 
(m) 

Radius 

Geoid 

(m) 

Difference 

(m) 


90N (pole) 

983.2231 

6356783.3 

6356796.3 

-13.0 


85 

983.1839 

6356944.9 

6356957.7 

-12.8 


80 

983.0676 

6357424.9 

6357437.1 

-12.2 


75 

982.8776 

6358208.9 

6358220.1 

-11.2 


70 

982.6195 

6359273.4 

6359283.2 

- 9.8 


65 

982.3011 

6360586.4 

6360594.6 

- 8.2 

N 

60 

981.9317 

6362108.6 

6362114.9 

- 6.3 


55 

981.5225 

6363793.9 

6363798.3 

- 4.4 

0 

50 

981.0857 

6365591.6 

6365593.9 

- 2.3 

R 

45 

980.6346 

6367447.2 

6367447.6 

- 0.4 

T 

40 

980.1829 

6369304.5 

6369303.1 

1.4 


35 

979.7444 

6371106.8 

6371104.0 

2.8 

H 

30 

979.3325 

6372799.4 

6372795.4 

4.0 


25 

978.9598 

6374330.3 

6374325.7 

4.6 


20 

978.6380 

6375652.7 

6375648.0 

4.7 


15 

978.3770 

6376726.0 

6376721.7 

4.3 


10 

978.1850 

6377517.1 

6377513.9 

3.2 


5 

978.0678 

6378001.8 

6378000.0 

1.8 

0 (equator) 

978.0300 

6378165.0 

6378165.0 

i 

0 


5 

978.0701 

6378001.8 

6378003.8 

- 2.0 


10 

978.1895 

6377517.1 

6377521.1 

- 4.0 


15 

978.3833 

6376726.0 

6376731.9 

- 5.9 


20 

978.6455 

6375652.7 

6375660.1 

- 7.4 


25 

978.9679 

6374330.3 

6374338.8 

- 8.5 

S 

30 

979.3404 

6372799.4 

6372808.3 

- 8.9 


35 

979.7514 

6371106.8 

6371115.4 

- 8.6 

0 

40 

980.1884 

6369304.5 

6369311.9 

- 7.4 

U 

45 

980.6379 

6367447.2 

6367452.8 

- 5.6 


50 

981.0862 

6365591.6 

6365594.7 

- 3.1 

T 

55 

981.5198 

6363793.9 

6363794.0 

- 0.1 

H 

60 

981.9257 

6362108.6 

6362105.4 

3.2 


65 

982.2919 

6360586.4 

6360579.9 

' 6.5 


70 

982.6074 

6359273.4 

6359263.7 

9.7 


75 

982.8629 

6358208.9 

6358196.4 

12.5 


80 

983.0509 

6357424.9 

6357410.3 

14.6 


85 

983.1660 

6356944.9 

6356929.0 

15.9 


90S (pole) 

983.2047 

6356783.3 

6356766.9 

16.4 


30 








With 


, . 1 . fa ^ 

2R J 
e 


the mean meridional oblateness becomes 


(40) 


f = 0.0033525944 = 1 : 298.276, 


(41) 


differing slightly from the Standard value of 1 : 298.30. 

The NASA Standard, here called Model 2, adopted by the Ad Hoc Com- 
mittee in 1963, has served and is serving its purpose well, which is to 
ensure accuracy and consistency between the various NASA agencies and 
contractors participating in the space program. While during 1961-63 
a real attempt was made to select a set of constants which might be 
termed the "best available at the time"^^, the chief qualification of 
the adopted set was standardization. Since the adoption of the NASA 
Standard for Apollo and other space programs, new and better determina- 
tions of the geographical (and astrodynamic) constants have been made. 
The Standard of 1963, still valid today, is therefore not based on more 
current information; however, with the present status of the Apollo pro- 
gram, it would probably be ill-advised to update the Standard at this 
time. 


The following models make use of the more recent determinations in 
various combinations. 


C. Model 3 (New * 13 , d 4 ) 


Similar to the preceding analysis, a model is assumed which takes 
the "pear" shape of the earth into account as well as the second even 
zonal harmonic, but not any dependence of the geopotential on longitude. 
The spherical harmonics coefficients are taken from Kozai^^ (equation 

( 13 )), 


J 

J 

J 


2 " 

1082.645 

3 " 

-2.546 

4 “ 

-1.649 


X 10- ® I 
X 10’® > , 
X 10"s 


(42) 
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where they have actually been determined together with ten additional 
coefficients, which are here assumed zero. This introduces a slight 
error of no practical concern for the present purpose. 

Gravity is again computed from equation (26), without the sectorial 
term, with J 2 , J 3 , J 4 as given in (42). 

According to the new J-coeff icients , the a- factor in the equation 
for the equatorial gravity acceleration must now be slightly adjusted; 
gg beccsnes now 


gg » 978.029 gal 

instead of the Standard value 978.030 gal. 

The reference ellipsoid (of revolution) postulated for this model is 
required to be consistent with the geopotential used. 

Using Krause' s theory, the A-coeff icients (equations (16) through 
(22)) in equation (15) are obtained. 


Aq = 0.9988804206 

Ai = 0 

As = 0.2236788408 x lO"^ \, 

A 3 = -0.2550678850 x 10"® 

A 4 = -0.3160571410 X 10‘S 


(43) 


and the radius of the assumed geoid becomes 

R 3 = 6371024.14 - 14266.61 P 2 (sin 0 ) + 16.268 P 3 (sin 0 ) + 20.16 P 4 (sin 0 ) 

(m) . (44) 

The results of this equation for the geoid surface are listed in 
column 4 of table 3. It can be seen that the polar radii are 

R „ = 6356794.0 m 
pN 

R _, = 6356761.4 m 
pS 
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TABLE 3 


Model 3 


Latitude 

(deg) 

Gravity 

(gal) 

Rad lus 
Ellipsoid 
(m) 

Radius 

Geold 

(m) 

Difference 

(Height of Ellipsoid) 
(n») 

90 

983.2239 

6356777.7 

6356794.0 

-16.3 

85 

983.1848 

6356939.4 

6356955.4 

-16.0 

80 

983.0685 

6357419.5 

6357434.8 

-15.3 

75 

982.8786 

6358203.7 

6358217.9 

-14.2 

70 

982.6205 

6359268.4 

6359281.2 

-12.8 

65 

982.3020 

6360581.8 

6360592.8 

-11.0 

N 60 

981.9326 

6362104.3 

6362113.3 

- 9.0 


981.5234 

6363790.1 

6363796.8 

- 6.7 

® 50 

981.0866 

6365588.3 

6365592.6 

- 4.3 

R 45 

980.6354 

6367444.4 

6367446.4 

- 2.0 


980.1835 

6369302.1 

6369302.0 

0.1 

^ 35 

979.7449 

6371105.0 

6371103.0 

2.0 

H 30 

979.3328 

6372797.9 

6372794.5 

3.4 

25 

978.9601 

6374329.3 

6374324.9 

4.4 

20 

978.6382 

6375652.0 

6375647.3 

4.7 

15 

978.3771 

6376725.6 

6376721.2 

4.4 

10 

978.1851 

6377516.9 

6377513.5 

3.4 

5 

978.0679 

6378001.7 

6377999.8 

1.9 

0 

978.0294 

6378165.0 

6378165.0 

0 

5 

978.0705 

6378001.7 

6378004.0 

- 2.3 

10 

978.1900 

6377516.9 

6377521.5 

- 4.6 

15 

978.3840 

6376725.6 

6376732.4 

- 6.8 

20 

978.6464 

6375652.0 

, 6375660.8 

- 8.8 

25 

978.9689 

6374329.3 

6374339.4 

-10.1 

S 30 

979.3416 

6372797.9 

6372808.8 

-10.9 

o 55 

979.7527 

6371105.0 

6371115.6 

-10.6 

° 40 

980.1896 

6369302.1 

6369311.8 

- 9.7 

U 45 

980.6389 

6367444.4 

6367452.2 

- 7.8 

T 50 

981.0871 

6365588.3 

6365593.4 

- 5.1 

^ 55 

981.5205 

6363790.1 

6363792.1 

- 2.0 

H 60 

981.9261 

6362104.3 

6362102.7 

1.6 

65 

982.2919 

6360581.8 

6360576.5 

5.3 

70 

982.6070 

6359268.4 

6359259.6 

8.8 

75 

982.8623 

6358203.7 

6358191.8 

11.9 

80 

983.0501 

6357419.5 

6357405.2 

14.3 

85 

983.1650 

6356939.4 

6356923.6 

15.8 

90 

983.2037 

6356777.7 

6356761.4 

16.3 
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so that the mean meridional flattening, equation (40), becomes 


f = 0.003353206 = 1 : 298.222. 


Comparison with the geoid of Model 2 shows that the new values for 
Jsi ^ 3 i J 4 cause the Model 3 geoid to be smaller at most latitudes (by 
2.3 m at the North Pole, 5.5 m at the South Pole, and by 1.2m and 0.6 m 
at 45N and 45S, respectively). 

For the reference ellipsoid, ^he equator radius imposed is the cur- 
rently accepted value due to Kaula, equation (35). With the above 
flattening, the polar radius of the ellipsoid then becomes 


R = R (1 - f) = 6356777.70 ra, 
p e' ' ’ 

and the ellipsoid shape can then be computed from equation (33). Results 
are listed in column 3 of table 3. The deviations of the ellipsoid from 
the geoid are given in column 5. The largest excursions are again at the 
poles, with 16.3 m height difference. 


D. Higher-Order Bodies with Equator Symmetry 
1. Model 4a (Triaxial Ellipsoid) 

In the preceding three ellipsoids, the longitude dependence of 
the geopotential has not been taken into account. For the present 
model, it is assumed that the geoid can be better approximated by a tri- 
axial ellipsoid. How well it agrees with the geoid is one of the objec- 
tives of this investigation. 

In a triaxial ellipsoid, the equator and all parallel sections 
exhibit a certain ellipticity, expressed -- in the geopotential -- by 
the sectorial coefficient J 2 2 * body is symmetrical with respect to 

the equator plane; hence, the odd-order terms in the spherical harmonics 
series of the geopotential (more exactly, of the radius vector of the 
surface) are zero. Thus, the earth is now assumed to be no longer 
symmetrical around the rotational axis. 
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The general formula of a triaxial ellipsoid is 



yf_ z2 


1 , 


( 45 ) 


where Reo Re^ are the two equatorial semi-axes. By introducing polar 
coordinates, Krause^"* has derived a general formula for the radius vector. 


cos^0 


^ + 4 ~ ^e ~ ^q) 

(1 - i 


ain^0 

(1 - f)2j 


lV2 


(46) 


where Rg is the mean equatorial radius, determined by the astrogeodetic 
and gravimetric methods of geodesy^* and currently assumed to be^ 


R « 6378165.0 m. 
e 


The symbol 0 is the geocentric latitude of the parallel, A is the longi- 
tude of the point on the surface, measured positive east of Greenwich, f 
is the meridional flattening of the ellipsoid, fg is the ellipticity of 
the equator (and all parallel sections), and Aq is the longitude of the 
apsidal line of the equatorial ellipse (i.e., the angle between the radius 
vector to the reference meridian on the equator and the semi-major axis 
of the equator ellipse, positive east of Greenwich), as well as the longi- 
tude of the associated parallel ellipses in both "hemispheres." Kaula® 
gives Aq = -14.5® ± 1.5®. In reference 1, he quotes -8® to -25®; how- 
ever, more recent data indicate an improved value of^*^ 

= -18® ± 3®. (47) 


Thus, the semi-major axis of the equator ellipse is rotated by 18® 
west of the Greenwich meridian. Consequently, the largest equator radius 
is located approximately 600 st. miles north of Ascension Island in the 
Mid-Atlantic, and - diametrically opposed - in the vicinity of the Solomon 
Islands in the Pacific, while the shortest radii are approximately 500 st. 
miles south and west of Ceylon in the Indian Ocean and about 1800 st. miles 
west of Peru in the Pacific Ocean. 
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By solving Krause's radius equation (equation (15)) twice for A = \, 
and A = Aq ± 90°, respectively, at the equator (0 = 0) and introducing 
the derived expressions for Reg and Rg^ in the definition of the flatten- 
ing (or oblateness). 


f 


e 



(48) 


one obtains a simple expression for the ellipticity of the equator, viz., 


f ~ 6*Ap 2 ~ Tj ds 2 ♦ 

The best current value for is given by Kaula®; 

Jo 2 = 1-8 ± 0.1 X 10"®. 


(49) 


(50) 


Using the most recent set of oblateness coefficients, equation (13), 
in the equation for the factor X, equation (22), one obtains 


X = 0.9981691858. 


With X and (equation (50)), equation (49) yields 


fg = 10.81980906 X 10"® = 1 : 92423. 


(51) 


(52) 


From equation (48) it is then found for the difference between the 
semi-major and semi-minor axis of the equator 


= Re - Re =89.0 m. (53) 

With the mean equatorial radius, Rg, according to equation (35), the 
longest radius of the equator becomes 

^e 

R + = 6378199.5 m 

e 2 
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and the shortest radius 


5 - ~ = 6378130.5 m 
e 2 

at the geographical locations described above. 

The sectorial harmonic is depicted in figure 1. 

The flattening f of the triaxial ellipsoid must be the same as for 
the "best" geoid (Model 5a and 5b) treated later, since the flattening 
equation, equation (58), contains only even zonal harmonics, in accord- 
ance with the symmetry behavior of the latter. With the even Jjj of 
equation (13), equation (58) yields 


f = 0.003353633015 = 1 ; 298.184. 


The polar radii then become 


R = (1 - f) R = 6378775.0 m. 

P e 

The surface of the triaxial ellipsoid has been computed for various 
latitudes and longitudes from equation (46). Results are given in 
table 4, column 2. 

The gravity of the Model 4a ellipsoid can be determined from equa- 
tion (26). All available even Jj, are used, equation (13), with the odd 
zonal coefficients set equal to zero according to the assumed symmetry 
with respect to the equator plane. With X as given in equation (51), 
the factor a becomes 0.9981683280, and the mean gravitation at the 
equator is then 


GM 

g = a — 978.0320 gal 

® R^ 
e 

differing by 2 mgal from Kaula*s Standard value. 
Results are presented in table 4, column 5. 
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TABLE 4 


Models 4a and 4b 


Longitude 

Radius of 

Rad ius of 

Height of 

Gravity of 

(deg 

)* 

Tri-ellipsoid 

Symm. Spheroid 

Spheroid 

Spheroid 

EOG 1 

WOG 

(m) 

(m) 

(m) 

(gal) 

LATITUDES 90N/90S (POLES) 

- 

6356775.0 

6356775.0 

0 

983.1815 

LATITUDES 75N/75S 


180 

6358203.0 

6358203.0 

0 

982.8385 

30 

150 

6358200.9 

6358200.9 

0 

982.8376 

60 

120 

6358199.0 

6358199.0 

0 

982.8367 

90 

90 

6358199.3 

6358199.3 

0 

982.8368 

120 

60 

6358201.4 

6358201.4 

0 

982.8378 

150 

30 

6358203.2 

6358203.2 

0 

982.8386 

LATITUDES 60N/60S 


180 

6362109.2 

6362114.2 

4.9 

981.9077 

30 

150 

6362101.4 

6362106.3 

4.9 

981.9041 

60 

120 

6362094.5 

6362099.3 

4.8 

981.9009 

90 

90 

6362095.4 

6362100.2 

4.8 

981.9013 

120 

60 

6362103.2 

6362108.1 

4.9 

981.9049 

150 

30 

6362110.1 

6362115.1 

4.9 

981.9081 

LATITUDES 45N/45S 

0 


6367456.9 

6367461.7 

4.7 

980.6259 

30 

150 

6367441.2 

6367446.0 

4.7 

980.6186 

60 

120 

6367427.4 

6367432.0 

4.6 

980.6122 

90 

90 

6367429.2 

6367433.8 

4. 6 

980.6131 

120 

60 

6367444.8 

6367449.6 

4.7 

980.6203 

150 

30 

6367458.7 

6367463.5 

4.8 

980.6267 

LATITUDES 30N/30S 

0 

B^l 

6372818.1 

6372819.3 

1.1 

979.3374 



6372794.6 

6372795.6 

1.0 

979.3263 

60 

120 

6372773.7 

6372774.7 

1.0 

979.3166 

90 

90 

6372776.4 

6372777.4 

1.0 

979.3179 

120 

60 

6372800.0 

6372801.0 

1.0 

979.3288 

150 

30 

6372820.8 

6372822.0 

1.2 

979.3384 


* EOG = East of Greenwich (Table 4 continued) 

WOG = West of Greenwich 
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TABLE 4 (Continued) 


Longitude 

Radius of 

Radius of 

Height of 

Gravity of 

(deg) 

Tri-Ellipsoid 

Symm. Spheroid 

Spheroid 

Spheroid 

EOG 

WOG 

(m) 

(m) 

(m) 

(gal) 

LATITUDES 15N/15S 

0 

180 

6376751.4 

6376749.3 

- 2.1 

978.3883 

30 

150 

6376722.0 

6376719.9 

- 2.1 

978,3747 

60 

120 

6376696.0 

6376693.9 

- 2.1 

978.3627 

90 

90 

6376699.4 

6376697.2 

- 2.2 

978.3643 

120 

60 

6376728.7 

6376726.6 

- 2.1 

978.3778 

150 

30 

6376754.8 

6376752.7 

- 2.1 

978.3898 

EQUATOR 

0 

180 

6378192.9 

6378192.9 

0 

978.0448 

30 

150 

6378161.4 

6378161.4 

0 

978.0303 

60 

120 

6378133.5 

6378133.5 

0 

978.0175 

90 

90 

6378137.1 

6378137.1 

0 

978.0192 

120 

60 

6378168.6 

6378168.6 

0 

978.0337 

150 

30 

6378196.5 

6378196.5 

0 

978.0465 


2. Model 4b (Symmetrical Geoid) 

With the oblateness coefficients up to J 14 readily available, 
it is of interest to investigate the deviation of the above triaxial 
ellipsoid from the surface of a quasi-ellipsoid, or rather "symmetrical 
geoid," derived from a series expansion of the known even zonal harmonics 
equivalent to the above employed gravity expansion. Using Krause's theory 
and solving equations (16) through (22), the coefficients in equation (15) 
are found to be 


Aq = +0.9988800610 
As = +0.2236730329 x lO"^ 
A 4 = -0.3180618021 X 10‘® 
As = -0.6471848752 x 10“® 
Aa = -0.2704952265 x 10 "® 
Aio= -0.5409904530 x lO’"^ 
Ai 2 = -0.3576547995 x lO'® 
Ai 4 = +0.1793283168 x 10"® 
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The surface equation, equation (15), then becomes 


R4 » 6371021.84 - 14266.235 - PsCsin 0) + 20.29 • P4(sin 0) 

- 4.13 • Ps(sin 0) + 1.73 . Pg(sin 0) + 0.345 • Pio(sin 0) 

+ 2.28 • Pi 2 (sin 0) - 1.14 • Pi 4 (sin 0) + 34.51 * cos^0 
» cos 2(A + 18 “) (m). (55) 


Equation (55) has been computed for some values of 0 and A; 
results are given in column 3 of table 4. The deviation of the sym- . 
metrical geoid from the true triaxial ellipsoid is listed in the fourth 
column of table 4. Agreement is seen to be extremely good. Height dif 
ferences are around 5 m in the higher latitudes. The largest height 
difference occurs around 55° North or South latitude, with over 6 m for 
all longitudes. 


E. "Best" Geoid Bodies 


1. Model 5a ("Best" Geoid with Sectorial Term) 

The most realistic model of the earth and its potential field 
must obviously make use of all available spherical harmonics coefficients, 
obtained from real-world observations. Using the theory developed 
recently by Krause^^> it is -- as shown above -- relatively simple to 
express the geoid surface, assumed an equipotential level, in an expan- 
sion of spherical harmonics and appropriate coefficients. This "true" 
geoid is assumed as Model 5a of this investigation. 

For the basic imposed constants (mean earth radius at equator, 
rate of the earth's rotation with respect to inertial space, and the 
gravitational parameter of the earth) the best currently available 
values^ are chosen: 

R = 6,378,165.0 m 
e 

GO = 0.729211585 x lO"'^ sec"^ 

GM = 398603,2 ± 3.0 km^ sec“^ 

© 
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and for the additional parameters 


Jo, 2 = 1.8 + 0.1 X 10"® 

\ = -18* ± 3*. 


Consistent with the above values, the mean equatorial gravita- 
tion is 


gg - 978.0320 gal (see Model 4a). 


Solving again equations (16) through (22) for all available J^, 
including odd orders, equation (13), the following An'coeff icients of 
equation (15) are obtained; 


"o 

Ai = 
Aa - 
A3 = 
A4 = 

As = 
As - 
Ay “ 
As “ 
Ag ® 

Aio ” 
All = 

Ai2 = 
Ai3 = 
Ai4 = 

with 

X = 
B2 = 


+0.9988800610 

0 

+0.2236730329, x lO"^ 
-0.2550669802 X lO'^ 
-0.3180618021 X 10-® 
-0.2103851762 x lO"® 
+0.6471848752 X 10"® 
-0.3336107793 x 10-® 
-0.2704952265 X lO"® 
-0.5309721113 X 10"^ ) . 

-0.5409904530 X 10"'^ 
+0.3025539200 X 10"® 
-0.3576547995 X 10"® 
-0.1142090956 X 10"® 
+0.1793283168 x 10"® 


0.9981691858 
0.3383437858 x lO"^ 


(56) 
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The equation for the surface of the geoid becomes then*; 


Rs = 6371021.844 - 14266.235 PsCsin j2!) + 16.268 P3(sin 0) 

+ 20.286 P4(sin 0) + 1.34 Ps(sin 0) - 4.127 PgCsin 0) 

+ 2.127 P 7 (sin 0) + 1.725 PgCsin 0) + 0.338 PgCsin 0) 

+ 0.345 Pio(sin 0) - 1.929 Pn(sin 0) + 2.28 PisCsin 0) 

+ 0.728 Pi 3 (sin 0) - 1.144 Pi 4 (sin 0)' + 34.505 cos^0 cos 2(A + 18®), 

(m) (57) 

where Pj^(sin 0) are the Legendre polynomials of argument sin p**. 

By writing the radius equation, equation (15), for the poles 
(0 = ±90°), one equation each results for the northern and the southern 
"hemispheres," which since 


f = 1 


R M + R o 
pN pS 


2R 


ead to an expression for the mean meridional oblateness of the geoid^^»^^ 




1 - (- 1 ) 


V 1.3.5. . . (2v - 1) 

2.4.6. . .2v 


v=l 


+ I Wg + B2(| Js - U)g - B 2 ) 

(58) 


} 


where X and are given in equation (56) and co^ = 3461.414 x 10"® sec"^ 
WithKozai's harmonics coefficients up to J 14 , equation (58) 

yields 


f = 0.003353633015 = 1 : 298.184, 


Some of the amplitudes (coefficients) of the geoid-equation in refer 
ence 14 have been found to be slightly incorrect. The values are 
corrected here. 

See Appendix A. 


42 



which is different from the currently accepted (Standard) value of 1960, 
viz. , 


f - 1 : 298.30, 

which is consistent with the old Jg (equation (37)) over a second-order 
Clairaut-type expression only, as mentioned earlier. 

The surface of the "best" geoid, equation (57), and its gravity, 
equation (26), have been computed and are tabulated in table 5 for various 
latitudes and longitudes. It can be seen that the geoid, in accordance 
with the fact that of all tesseral harmonics only the first sectorial 
was included, is symmetrical around a plane defined by the major axis 
(at * -18*)~'of the equator ellipse and the axis of rotation. If 
higher tesserals were taken into account, this last symmetry would also 
disappear. 

The sectorial harmonic at the equator is shown in figure 1. The 
curve depicts (greatly exaggerated) the height of the "best" geoid over 
the rotationally symmetric geoid of Model 4b around the equator. At 
latitudes 45M and 45S, the amplitudes of the harmonic, for equal meridians, 
are half their value at the equator. 

Because of the sectorial term depicted, the longest equator radius 
(also the longest radius vector to all latitudes) is at 162°E and at 18°W 
of Greenwich, with 6,378,199.5 m. The shortest radial vector is at 72®E 
and 108 °W of the prime meridian, with 6,378,130.5 m. The difference 
amounts to 69.0 m. 

At 45 degrees North and South parallel, the difference between 
the semi-major and semi-minor axis, as stated above, must be one-half of 
this value. It is 34.5 m. 

At the poles, it is seen that the radial vector at the North Pole 
is longer than at the South Pole. The difference in height is 37.8 m, 
which accounts for the "pear" shape of the earth, suspected already in 
1959 by O'Keefe®. The gravity at the North Pole is higher by 13.11 mgal 
than at the South Pole. Around the equator belt, the gravity varies by 
as much as 31.74 mgal (between the longitudes 72® and 162®E), The mean 
equatorial gravity must be the same as for the triaxial ellipsoid. Model 4a; 
i.e. , 


g^ = 978.0320 gal. 
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TABLE 5 
Geold 


Longitude* 
EOG 1 WOG 

Radius of 
Geo id 

(m) 

Gravity of 
Geo id 

(gal) 

Longitude* 
EOG 1 WOG 

Radius of 
Geoid 
(m) 

Gravity of 
Geoid 
(gal) 




UTITUDE 90N 

(North Pole) 



- 

6356793.9 

983.18807 


LATITUDE 75N 

0 

180 

6358217.6 

982.84373 


90 

6358213.8 

982.84200 

30 

150 

6358215.4 

982.84275 


60 

6358215.9 

982.84298 


120 

6358213.6 



982.84189 

150 

30 

6358217.8 

982.84384 

LATITUDE 60N 

0 

180 

6362118.1 

981.90820 



6362104.2 

981.90178 

30 

150 

6362110.3 

981.90457 


60 

6362112.1 

981.90540 

60 

120 

6362103.3 

981.90136 

150 

30 

6362119.0 

981.90861 

LATITUDE 45N 

0 

iil 

6367458.7 

980.62507 

WM 

90 

6367430.8 

980.61223 


wm 

6367443.0 

980.61782 

mm 

60 

6367446.6 

980.61948 


120 

6367429.0 

980.61140 


30 

6367460.5 

980.62590 

LATITUDE 30N 

0 


6372812.7 

979.33540 

90 

90 

6372770.8 

979.31614 

30 

■BW 

6372789.0 

979.32452 

120 

60 

6372794.4 

979.32701 


120 

6372768.1 

979.31489 

150 

30 

6372815.4 

979.33664 

LATITUDE 15N 

0 

180 

6376743.6 

978.38634 

90 

90 

6376691.6 

978.36238 

30 

150 

6376714.2 

978.37281 

120 

60 

6376721.0 

978.37591 

60 

120 

6376688.2 

978.36083 

150 

30 

6376747.0 

978.38789 


(table 5 continued) 

'k 

EOG = East of Greenwich (degrees) 

WOG = West of Greenwich (degrees) 
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TABLE 5 (Continued) 


Longitude 
EOG I WOG 


Radius of 

Gravity of 

Longitude 

Radius of 

Geold 

Geold 

EOG i WOG 

Geold 

(m) 

(gal) 

(m) 


LATITUDE 0 (EQUATOR) 



6378192.9 

978 04484 

6378161.4 

978.03034 

6378133.5 

978.01750 



6378137.1 

978.01916 

6378168.6 

978.03366 

6378196.5 

978.04650 


LATITUDE 15S 



6376755.0 

978.39020 

6376725.6 

978.37667 

6376699.6 

978.36469 



6376702.9 

978.36623 

6376732.3 

978.37976 

6376758.4 

978.39174 


LATITUDE 30S 


0 

180 

6372825.9 

979.33889 

30 

150 

6372802.2 

979.32801 

60 

120 

6372781.3 

979.31838 


6372784.0 

6372807.7 

6372828.6 

979.31963 

979.33050 

979.34013 




LATITUDE 45S 



6367464.7 

980.62674 

6367448.9 

980.61949 

6367435.0 

980.61307 



6367436.8 

980.61390 

6367452.5 

980.62115 

6367466.5 

980.62757 


LATITUDE 60S 



981.9 

981.9 

981.9 



M 

90 

6362096.3 

981.90080 

60 

6362104.1 

981.90443 

150 

30 

6362111.1 

981.90764 


LATITUDE 75S 



6358188.5 
6358186.3 

6358184.5 

982.83336 

982.83239 

982.83153 

LATITUDE 9 

1 6356756.1 

983.17496 



6358184.7 

982.83164 

6358186.8 

982.83261 

6358188.7 

982.83347 
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2. Model 5b (”Best" Geoid without Sectorial Term) 

For purposes of comparison with models fea*curing circular equa- 
tors, it is also of interest to determine the shape of the **best” geoid, 
with all available zonal harmonics coefficients, based on a circular 
equator. Model 5b is therefore obtained by solving equation (57) with- 
out the last (sectorial) term. Since the sectorial harmonic at the equa- 
tor has already been determined in figure 1 , the desired circular-equator 
shape is obtained by subtracting it from the radius vectors of Model 5a. 


IV. COMPARISON OF MODELS 


In accordance with the different assumptions used, the various geo- 
graphical theories to be compared here define geometrical bodies which 
differ from each other and from the geoid. Also, they lead to different 
gravity values for identical locations on the body. Primary attention 
is directed towards the deviations in geometries. 

Figure 6 shows the radius vector as a function of latitude of Models 
1 , 2 and 3, referenced to the same body, i.e., the ’’best” geoid without 
the sectorial term (Model 5b), In all three cases, the bodies are not 
ellipsoids, but spheroids of revolution symmetrical around the equator 
plane, since they are based on superimposed spherical harmonics. It is 
seen that the widest deviation from the geoid is exhibited by Model 1, 
which takes into account solely the second order zonal harmonic, At 

the North Pole, the difference is -32,6 m; at the South Pole, +5.20 m. 
Since the spheroid of Model 1, as shown on table 1, differs from the 
Model 1 ellipsoid (flattening 1/298.254) by as much as -18,7 m at the 
North and South Pole, the deviation of the Model 1 ellipsoid from the 
reference geoid is less at the North Pole, amounting to -13,9 m, and 
higher (23.9 m) at the South Pole. 

The spheroids of Models 2 and 3 are based on the old and the new 
set of coefficients J^, J 3 , J 4 , respectively, as given by equations (39) 
and (44). In general, they agree very well with the reference geoid, 
exhibiting deviations of less than 5 m. At the poles, the deviation of 
the Model 2 spheroid from the reference line as well as from the improved 
Model 3 spheroid is larger, primarily due to the changes in and J 3 . 

As in the case of Model 1, it is of interest to determine .the devia- 
tion of the Model 2 spheroid from its associated ellipsoid, as given in 
table 2. The height of the spheroid based on the Standard values of J 2 , 
J 3 , J 4 over the ellipsoid with consistent flattening 1/298.30 is shown 
in figure 7. Also plotted is the height of the "best” geoid without the 
J 2 2 "t®™ (Model 5b). The agreement between the two Model-2 bodies is 
good; the largest difference is at the South Pole, with 16.4 m. As is 
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South Polo 


6. 6. HEIGHT OF MODELS i, 2 AND 3 (SPHEROIDS) 
OVER ‘BEST* 6E0ID WITHOUT SECTORIAL TERM 


47 






North Pole 



South Pole 

FIG. 7. HEIGHT OF MODEL 2-SPHEROID AND OF "BEST" GEOID 
WITHOUT SECTORIAL TERM, OVER STANDARD ELLIPSOID ('1/298.3) 
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to be expected, the spheroid agrees better with the "best" geoid than the 
ellipsoid. The influence of particularly the odd zonal harmonics becomes 
manifest at the South Pole, where the Standard Ellipsoid, as used pre- 
sently by NASA, deviates by as much as 27.2 m from the Model 5b-geoid. 

A still better approximation to the "best" geoid with circular equa- 
tor can be given with the spheroid of Model 3, as depicted in figure 8 . 
Here, the spheroid height is referenced to the associated ellipsoid of 
flattening 1/298.222, derived from the set of improved J2, J3, and J 4 . 

The deviations of the spheroid from the ellipsoid are numerically very 
similar to those between the Model 2 bodies, while the agreement of the 
ellipsoid with the circular-equator geoid has been improved considerably, 
except for the latitudes around 20N and SOS, where the gap is still almost 
10 m. The agreement of the spheroid with be "best" geoid is generally 
better than 5 m. 

The height of the "best" geoid (Model 5a) has been computed from the 
data of tables 4 and 5, and is presented in figure 9 (greatly exaggerated), 
referenced to the triaxial ellipsoid of Model 4a. Both bodies are based 
on expansions of spherical harmonics up to and including 3 ^ 4 , with the 
latest available values for the coefficients. The sectional cut is along 
the prime meridian. The "pear" shape of the earth is clearly discernible. 
Since the mean equatorial radius is the same for both bodies (as well as 
the sectorial harmonic), they share the same equator. At the poles, how- 
ever, deviations of the geoid are observable, amounting to 18.87 m at the 
North Pole, and to -18.87 m at the South Pole. In the intermediate lati- 
tudes, the deviations are less, with -8.5 m at 20 N and 7.7 m at 30S. The 
height differences are also plotted in figure 10 , on a rectilinear grid, 
for a meridional section at Greenwich longitude. 

According to the different harmonics coefficients used in the poten- 
tial equation, the gravity of the analyzed models differs in each case. 
Using the gravitation of the "best" geoid (Model 5a) as tabulated in 
table 5 , as common reference, the gravities of each model were computed 
for a meridional section along the prime meridian and plotted in figure 11 . 
As is to be expected, the gravity of the triaxial ellipsoid (Model 4a), 
listed in table 4, agrees best (generally by better than ±3 mgal) with 
the "best" geoid. The other models, with gravity values determined from 
differing values and numbers of oblateness coefficients, show generally 
identical behavior, with significant deviations (up to 35 mgal) at the 
poles and equator, in accordance with the recent findings of a more pro- 
nounced influence of the odd zonal harmonics, particularly J 3 . 

For Model 1, the exclusive use of Jg apparently compensates somewhat 
around the North Pole for the deviations in Models 2 and 3, which are 
obviously caused by the third order zonal harmonic, the only additional 
additive amplitude at these latitudes besides Jg. At the lower latitudes. 
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FIG. 8 HEIGHTS OF MODEL 3- SPHEROID AND 'BEST* GEOID WITHOUT SECTORIAL TERM 
OVER ELLIPSOID OF FLATTENING 1/298.222 (MODEL 3) 
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North Pole 



FIG. 10. HEIGHT OF 'BEST* GEOID (MODEL 5) 

OVER TRIAXIAL ELLIPSOID OF FLATTENING 1/298.184, 

ALONG PRIME MERIDIAN 
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Model 3 — with the improved set of coefficients — agrees best of all 
three models with the reference, as is to be expected. However, the 
improvement is only very minor (< 1 mgal), leading to the observation 
that the effects of improving the oblateness coefficients (rather than 
increasing their number beyond J 4 ) are apparently revealed more signi- 
ficantly in the geometries than in the gravities of the bodies used. 


V. CONCLUSIONS 


1. A new earth figure model, or geoid, has been developed^^* and 
computed here, consistent with the most recent oblateness coefficients 
available. 

2. The currently used NASA Standard ellipsoid differs from this 
"best" geoid by -10.6 m at the North Pole, 9.6 m at 20°N parallel, -4.6 m 
at 35®S parallel, and as much as 27.2 m at the South Pole. If a spheroid 
based on a series expansion of spherical harmonics with the same oblate- 
ness coefficients as in the gravity equation were used instead, as pro- 
vided by Krause's theory, these deviations could be reduced to 2.4 m, 

4.9 m, 4.0 m, and 10.8 m, respectively. This surface formula would be 
more consistent with the gravity formulation. 

3. The currently used NASA Standard oblateness coefficients J 2 , 

J 3 , J 4 result in gravity values, which differ from the gravity based on 
recent satellite measurements by as much as 35 mgal at the North Pole, 

-15.7 mgal at the equator, and 30 mgal at the South Pole (along the 
Meridian), Replacing the three J-coeff icients by their updated values 
does not lead to significant improvements. 

4. Conversely, updating the oblateness coefficients Jg, J 3 , J 4 in 
the NASA Standard geometry model (and using the series expansion, as noted 
above) could lead to significant improvements in height deviations (see 
figure 6 ), amounting to about 100 percent at the North Pole, 50 percent 

at the South Pole, and 30-50 percent at most of the intermediate latitudes. 

5. For certain applications (e.g., guidance of short duration first- 
stage boosters), use of geographical and gravitational models based only 
on J 2 may be permissible. Because of the shape of the Jg-harmonic (fig- 
ure 2), it is advisable to use the ellipsoid (flattening 1/298.254), not 
the spheroid, as geometrical reference. Its height deviations from the 
"best" geoid (circular equator) are less than 10 m in the belt between 

65 ®N and 65®S latitudes. The gravity of this simple model differs by 
maximally 30 mgal, at the South Pole (-18.5 mgal at equator, 16.8 mgal 
at 90®N) from the sophisticated gravity model. The deviation is zero at 
about 24®N and 30®S parallels intersecting the Prime. 
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6. The latest available satellite measurements result In a **best** 
geold which Is pronouncedly "pear” shaped, with the north polar radius 
being 37.8 m longer than the south polar radius vector. The equator is 
slightly elliptic, with a difference of 69.0 m between the semi-major and 
semi-minor axes. The major axis of the equator ellipse is rotated West 
against the Prime Meridian. The rotation angle (longitude) is currently 
thought to be -18® ± 3®. As shown by Krause, the "best" geold can be 
described by a series expansion of spherical harmonics of the same degree 
as in the geopotential, the (amplitude) coefficients of which are directly 
related to the satellite-obtained oblateness coefficients. 

7. As an approximation to the "best" geoid, an exact ellipsoid can 
be defined, as is customary in geodesy. The ellipsoid is a triaxial 
ellipsoid with the same equator ellipse as the geoid. The largest devia- 
tions of this body from the geoid are 18.87 m at the poles. In a band 
bounded by the parallels 63®N and 70®S, the height differences amount to 
less than 10 m. If for the geopotential of the triaxial ellipsoid only 
even J's are used, the gravity deviates from the geopotential of the geold 
generally by less than 7 mgal and between 70®N and 70®S by less than 5 mgal. 
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APPENDIX A 


Computation of Legendre Polynomials 


The Legendre spherical functions P^C^^) polynomials of the type 


P 


n 



3*5. 


(2n-3H2n-l) 

nl 


n_ n(n-l) n -2 
2(2n-l) ^ 


nCn-l)(n-2)(n-3) 

2-4(2n-l)(2n-3) 


n-4 

X 


• • • 


n = 0, 1, 2 

They are solutions of the differential equation 
(1 - x2) 0 - 2x g + n(n + 1) y = 0. 


For the present analyses, the Legendre functions were computed from 
the recurrence relation 


(n + 1) - (2n + 1) X P^(x) + nP^_^(x) = 0. 


where 


X = sin 


and starting with 


PgCsin 0 ) = 1 

Pj(sin 0) - sin0. 


The most extensive tabulation of Pj^(cos 0) , for n = 0(1)80 and 
0 = 0“(1®)180‘’, has been published in reference 20. Since 


cos 0 = sin (90® ± 0) , 


the spherical functions can also be taken from this reference. 


pRECEomci umR noT nwss>. 
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APPENDIX B 


Alternate Expression for the Gravity Acceleration 


Krause, in references 13 and 14, gives a form of the gravity. 



(n + 1) Cj^P^(sin 0) + 30^^ cos 2(A “ A^j) 


( 59 ) 


obtained from equation (26) by several transfoixoations^^. The coef- 
ficients Cn are related to the oblateness coefficients and can be 
determined from them by relations similar to the equations (16) through 
(22) for the coefficients An, involving some of the latter. 


C 

o X 


1 + 2(1 - Aq) ■ f “e ■ J BaClSJa + Wg - 3Ba) 


C T ""0 


Cp = “ 

^ X 


J2 - 3 As - 9 Wg - 


(3J2 + Wg)^ + -^ B2(18J2 + - 3B2) 


C4 = T7 


J 4 - y A 4 + YYj (3J2 + ^ 175^ B2(18J2 + ~ 3B2)J 


2A, 


X 


J„ - 


K K + 1 


for K ^ 0,2,4 


O2. 2 ~ ^ J2. 2 


For the "best" geoid (Model 5a) with Kozai's set of Jn up to J14, 
the C-coef f icients become 


Cq = 1.001760646 

Cl = 0 

C2 = -0.1172324264 x lO’^ 


pREes»<j m€S pfOT nmeo. 
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C 3 = -0.1272995717 X lO"^ 
C 4 =+0.2636592115x10-® 
C 5 = -0.1401281566 x 10"® 
Cs = +0.4619357542 x 10"® 
Cy = -0.2500551096 X 10'® 
Cq = -0.2102749240 X 10‘® 
C 9 = -0.4245829102 X lO""^ 
Cio = -0.4424481398 X lO”^ 
Cii = +0.2520357773 x 10“® 
C 12 = -0.3025300610 X 10"® 
Ci 3 = -0.9786358495 x lO’’^ 
Ci 4 = +0.1553740187 x 10"®. 


With these coefficients, equation (59) then becomes*, using 
gg = 978.0320, 

g = 979753.968 + 3439.712 PgCsin 0) + 4.980 P3(sin 0) - 12,893 P4(sin 0) 
+ 0.822 PsCsin 0) - 3.163 Ps(sin 0) + 1.956 P^Csin 0) 

+ 1.851 Pa(sin 0) + 0.415 PgCsin 0) + 0.476 Pio(sin 0) 

- 2.958 Pii(sin 0) + 3.846 P^gCsin 0) + 1.340 Pi3(sin 0) 

- 2.279 Pi 4 (sin 0) + 15.873 cos^0 cos 2(A+ 18°) (mgal), (60) 


With the Legendre polynomials obtained as shown in Appendix A, the 
gravitation can be computed for the "best" geoid from equation (60) as 
an alternative to the gravity equation (26). 




ii' 


Some of the amplitudes (coefficients) of the gravity equation in 
reference 14 have been found to be slightly incorrect. The values are 
corrected here. 
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